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Abstract. Dynamics of inelastic gases are studied within the framework of random 
collision processes. The corresponding Boltzmann equation with uniform collision 
rates is solved analytically for gases, impurities, and mixtures. Generally, the energy 
dissipation leads to a significant departure from the elastic case. Specifically, the 
velocity distributions have overpopulated high energy tails and different velocity 
components are correlated. In the freely cooling case, the velocity distribution de- 
velops an algebraic high-energy tail, with an exponent that depends sensitively on 
the dimension and the degree of dissipation. Moments of the velocity distribution 
exhibit multiscaling asymptotic behavior, and the autocorrelation function decays 
algebraically with time. In the forced case, the steady state velocity distribution 
decays exponentially at large velocities. An impurity immersed in a uniform inelas- 
tic gas may or may not mimic the behavior of the background, and the departure 
from the background behavior is characterized by a series of phase transitions. 

1 Introduction 

Granular gases consist of weakly deformable macroscopic particles that in- 
teract via contact interactions and dissipate energy during collisions [1]. In 
typical experimental situations the system is constantly supplied with en- 
ergy to counter the loss occurring during inelastic collisions [2]. The most 
ubiquitous features of granular gases are their tendency to form clusters and 
their anomalous velocity statistics [3-6]. In turn, these characteristics lead to 
many interesting collective phenomena including density inhomogeneities [7] , 
shocks [8,9], size segregation [10], pattern formation [11], and phase transi- 
tions [12], to name a few. 

Inelastic hard spheres provide a useful theoretical, computational, and 
experimental framework for studying granular gases [13-16]. Inelastic gases 
pose new theoretical challenges as their dynamics are governed by dissipative 
collisions and strong velocity correlations. Fundamental equilibrium charac- 
teristics such as ergodicity and equipartition of energy are typically violated 
by this nonequilibrium gas system. As a result, inelastic gases exhibit fi- 
nite time singularities [17-19], chaotic behavior [20], breakdown of molecular 
chaos [21-25], and anomalous velocity statistics [26-30]. 

The significance of the energy dissipation is nicely demonstrated by con- 
sidering a mean-field version of the hard sphere system, namely a random 



2 Ben-Nairn and Krapivsky 



collision process. In this process, randomly chosen pairs of particles undergo 
inelastic collisions with a random impact direction. This process, often called 
the Maxwell model, is described by a Boltzmann equation with a velocity 
independent collision rate. As this model is random in all respects, it is an- 
alytically tractable. In fact, in classical kinetic theory of gases, the Maxwell 
model precedes the Boltzmann equation [31]. Historically, it played an im- 
portant role in the development of kinetic theory [32-36] , and it remains the 
subject of current research [37-39]. Recently, it has been realized that the 
Maxwell model is analytically tractable even for inelastic collisions [40-55]. 
In this review, we detail dynamics of uniform inelastic gases, isolated impu- 
rities in uniform gases, and binary mixtures. 

First, we consider the onc-dimcnsional gas where explicit analytic solu- 
tions are possible (Sec. 2). In the freely cooling case, the velocity distribution 
approaches a universal (dissipation independent) scaling form with an alge- 
braic high-energy tail. Furthermore, the moments of the velocity distribution 
exhibit a multiscaling asymptotic behavior, and the velocity autocorrelation 
function decays algebraically with a non- universal exponent. In the uniformly 
heated case, the system approaches a steady state, and the velocity distribu- 
tion has an exponential high-energy tail. 

Many of these features extend to higher dimensions (Sec. 3). In the freely 
cooling case, however, the scaling function underlying the velocity distribu- 
tion is no longer universal — its shape and its extremal behavior depend on 
the spatial dimension and the degree of dissipation. We obtain explicit ex- 
pressions for low order moments and show the development of correlations be- 
tween different velocity components. Such correlations are more pronounced 
in the freely cooling case, compared with the forced case. 

When an impurity is immersed in a homogeneous cooling fluid, various 
scenarios are possible (Sec. 4). If the impurity mass is lighter than a critical 
mass, its velocity statistics are completely governed by the background. As 
the impurity mass increases, velocity moments of higher order may diverge 
asymptotically compared with the fluid. A series of critical masses govern 
these phase transitions. Ultimately, the temperature ratio may diverge, and 
the impurity becomes so energetic that it effectively decouples from the fluid. 
In this case, its velocity distribution is distinct from the fluid. Although single- 
time velocity statistics of the impurity may be governed by the fluid, two-time 
statistics such as the autocorrelation function are always different. 

For mixtures (Sec. 5), all components have the same asymptotic temper- 
ature decay, and a high-energy tail as in the uniform gas case. This behavior 
is shown in detail for one-dimensional freely cooling binary mixtures. 

We finally describe a lattice generalization of the Maxwell model where 
particles are placed on a lattice and only nearest neighbors interact via inelas- 
tic collisions (Sec. 6). We show the development of spatial velocity correlations 
with a diffusively growing correlation length. These correlations significantly 
reduce the temperature cooling rate. 
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2 Uniform Gases: One Dimension 
2.1 The freely cooling case 

Consider an infinite system of identical particles that interact via instan- 
taneous binary collisions. When a pair of particles with velocities (1x1,1*2) 
collide, their post collision velocities (vi,V2) are given by the collision rule 

"1,2 = wi,2 T (1 -p)g, (1) 

with g = mi — ii2 the relative velocity and p the collision parameter. The total 
momentum is always conserved, u\+u 2 = vi+v 2 - Since the system is invariant 
under the Galilean transformation v — > v — v , we set the average velocity to 
zero, without loss of generality. In each collision, the relative velocity between 
the colliding particles is reduced by the restitution coefficient r = 1 — 2p. In an 
inelastic collision, there is an energy loss, AE = — p(l — p)g 2 (particle mass 
is set to unity); energy loss is maximal for completely inelastic collisions 
(p = 1/2), and it vanishes for elastic collisions (p = 0). 

In a random collision process, the collision rate is independent of the 
particles' velocities. Let P(v,t) be the normalized velocity distribution at 
time t. It evolves according to the Boltzmann equation 

dP(v,t) =R r duiP(uijt) f dU2 p( U2it) (2 ) 



dt 

x {5[v-U! + (l- p)g] - S(v - u 2 )} . 

This Boltzmann equation, with a velocity independent collision kernel, is 
termed the Maxwell model in kinetic theory [35]. The overall collision rate 
K = with the granular temperature T = J dv v 2 P(v, t), is chosen to 
represent the typical particle velocity. The quadratic integrand in Eq. (2) 
reflects the random and binary nature of the collision process, while the gain 
and loss terms reflect the collision rules (1). One can verify that the total 
number of particles, J dvP(v, t), and the total momentum, J dv vP(v, t), are 
conserved; hereinafter we set / dvP(v,t) = 1 and / dvvP(v,t) = 0. 

For inelastic hard spheres, where the collision rate equals the relative 
velocity, the rate equation governing the temperature is part of an infinite 
hierarchy of equations. In contrast, in the Maxwell model, the temperature 
obeys a closed rate equation 

dT 

_ = _ XT 3/2 with \ = 2p(l-p). (3) 
The temperature decays algebraically with time according to Haff's law [56] 



(4) 
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with the time scale t^ 1 = p(l — p)v / To- Eventually, all particles acquire the 
average velocity: P(v, oo) = S(v). Our goal is to characterize how the velocity 
distribution approaches this terminal state asymptotically. Since the overall 
collision rate is uniform, it is useful to characterize time by the collision 
counter r = J dt'K(t'), equal to the average number of collisions experienced 
by a particle, 

t = j ln(l + t/t ). (5) 

In terms of the number of collisions, T(r) = To e~ Ar . 

The convolution structure of the Boltzmann equation suggests to apply 
the Fourier transform F(k,r) = J dv e lkv P{v,t). This quantity evolves ac- 
cording to [40] 

|- F(k, t) + F(k, t) = F(k - pk, t) F(pk, r). (6) 

This closed equation is both nonlinear and nonlocal, yet it is analytically 
tractable. Assuming that the velocity distribution approaches its final state 
in a self-similar fashion, we seek a scaling solution 

P(v,t) =T- 1/2 T(w) with w = vT- 1/2 . (7) 

The scaling form corresponding to the Fourier transform is F(k, r) = f(z) 
with the variable z = I&IT 1 / 2 . Substituting the scaling ansatz into Eq. (6) 
and using dT/dr = —AT, the Fourier scaling function satisfies 

-p(l - p)f'(z) + f(z) = f(z - pz)f{pz). (8) 

This equation is supplemented by the small-z behavior f(z) = 1 — \z 2 , which 
is dictated by the small wave number behavior, F(k) = 1 — \k 2 T . Subject to 
these conditions, the (unique) solution is [42] 

f(z) = (l + z)e- z . (9) 

The scaled velocity distribution is obtained by performing an inverse Fourier 
transform 

V(w) = - 1 (10) 

7T (1 + W ) 

Remarkably, the scaled velocity distribution is independent of the collision 
parameter p. Another important feature is the algebraic decay: V(w) <~ u>~ 4 
for w 3> 1. This behavior should be compared with the exponential high 
energy tails obtained for the traditional Boltzmann equation [26,28]. The 
enhancement in the likelihood of finding energetic particles is due to the 
effective reduction in their collision rate. 

Typically, in kinetic theory of molecular gases, the velocity distributions 
have sharp tails such that all moments of the distribution are finite 1 . Sonine 

1 Algebraic high-energy tails may also characterize nonequilibrium states of elastic 
gases; for uniform shear flows of two-dimensional Maxwell molecules, this has 
been recently shown by Acedo, Santos, and Bobylev [38]. 
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expansions, where Maxwellian distribution are systematically modified by 
polynomials of increasing orders, are widely used to analyze the Boltzmann 
equation. The Maxwell model can be conveniently utilized to examine the 
applicability of this approach to situations where the velocity distributions 
have ovcrpopulated high energy tails. 

The scaling function does not characterize all features of the asymptotic 
time dependent behavior. Higher than third moments of the scaling function 
diverge, and since moments must be finite at all times, a direct calculation 
is necessary. Moments of the velocity distribution, M n (t) = J dvv n P(v,t), 
obey a closed set of equations 

— M n + a n M n =Y.[ W 1 - p)"- m M m M„_ m , (11) 

with a n (p) = 1 — p n — (1 — p) n - These equations are solved recursively using 
Mo = 1 and Mi = 0. Assuming that all moments are initially finite we find 
that to leading order, the moments asymptotically decay as 

M n ~ e~ a " T - t- 2a "? a2 . (12) 

Indeed, the sum term in (11) is asymptotically negligible because the coef- 
ficients satisfy the inequality a n < a m + a„_ m for all 1 < m < n — 1 [40]. 
Thus, the temperature does not characterize higher order moments. Instead, 
M n ~ M£ n with a nonlinear multiscaling spectrum a n = a n /a2, rather than 
the linear spectrum a n = n/2 that characterizes ordinary scaling behavior. 
The indices a„ increase monotonically with p, so the stronger the dissipation, 
the more pronounced the multiscaling asymptotic behavior. 

The scaling behavior characterizes only velocities of the scale of the typi- 
cal velocity. For sufficiently large velocities, far outside the scaling region, the 
gain term in the Boltzmann equation is negligible, and ( + l)P(v, t) = 0. 
Consequently, a generic exponential decay, P(v,t) ~ Pq(v) exp(— t) charac- 
terizes such large velocities. Similarly, sufficiently large moments decay ac- 
cording to M n <~ e~ T <~ i _1 /[p( 1_ p)l for n — > oo. 

In several studies, inelasticity is treated as a small perturbation [29,30,57]. 
To first order in p, Eq. (6) reads (J^ + pkJ^) F(k,r) = 0. The solution to 
this equation, F(k,r) = Fo(ke~ pT ), remembers the initial conditions forever, 
in contradiction with the exact asymptotic behavior (9) . This example raises 
questions concerning the validity of perturbation analysis in the vicinity of 
p = 0. The p — * limit is singular and indeed, in Eq. (8) the small parameter 
p multiplies the highest derivative. 

The Fourier transform can be expressed as a series expansion. The linear 
term can be eliminated from Eq. (6) by making the transformation F — > e T F 
and r — ► 1 — e~ r . Then a formal Taylor expansion solution can be found: 

F(k,T) = e-J2 [ - r^F^k). (13) 
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The expansion functions F n (k) are obtained recursively: F (k) = F(k, r = 0), 
and generally F n+1 (k) = YJL=o (Z) F m(k-pk)F n ^ m (pk). For instance, when 
the initial distribution is Maxwellian, the expansion functions, and hence the 
velocity distribution itself consist of sums of Maxwellians. Similarly, start- 
ing from a stretched exponential, all expansion functions consist of sums of 
stretched exponentials. In general, the expansion functions are products of 
F 's with stretched arguments. This implies that starting from a compact 
initial distribution P(v,0), the velocity distribution P(v,t) develops a set of 
singularities. For instance, a distribution with support in [— vq,vq} becomes 
non-analytic at an infinite set of points vi^ m = ±p l (l — p) m vo- 

Thus far, we characterized velocity statistics at a specific point in time. 
The autocorrelation function 

A(t',t) = v(t')v(t) (14) 

with the overlinc denoting an average over all particles, quantifies (two-point) 
temporal correlations in the velocity of a tagged particle. We have 

^-A{t',t) = v{f) [du(r)/dr] = -(1 - p) u(r') [v(t) - u(t)} 

= -(1- P )A(t',t). (15) 

Therefore, the autocorrelation decays exponentially with the collision num- 
ber. In terms of the original time variable: 

A(t', t) = A (1 + t'/t ) 1/p ~ 2 (1 + t/to)~ 1/p , (16) 

with Ac = Tq. In particular, A(t) = A(0, t) oc t _1 / p , so memory of the 
initial conditions decays algebraically with time. The autocorrelation function 
decays faster than the temperature, A(t) < T(t) (the two functions coincide 
in the completely inelastic case, p — 1/2). Generally, the autocorrelation is a 
function of the waiting time t' and the observation time t, and not simply of 
their difference, t — t'. This history dependence ("aging") merely reflects the 
fact that the collision rate keeps changing with time. 

The spread in the position of a tagged particle, A 2 (t) = (\x(t) — x(0)| 2 ), is 

obtained from the autocorrelation function using A 2 = 2 f* dt' J* dt"A(t", t'). 
Substituting the autocorrelation function (16), we find that asymptotically, 
the spread grows logarithmically with time 

A oc Vlai. (17) 

This behavior reflects the t~ l decay of the overall velocity scale [58,59]. 

2.2 The forced case 

In experimental situations, granular ensembles are constantly supplied with 
energy, typically through the boundaries, to counter the dissipation occurring 
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during collisions 2 . Theoretically, it is convenient to consider uniformly heated 
systems, as discussed by Williams and MacKintosh [14] in the case of one- 
dimensional inelastic hard rods and then extended to higher dimensions (see 
[28-30] and references therein) . We therefore study a system where in addition 
to changes due to collisions, velocities also change due to an external forcing: 
-fo | heat = £j- We assume standard uncorrelated white noise: = and 

Such white noise forcing amounts to diffusion with a 'diffusion' coefficient 
D in velocity space. Therefore, the Boltzmann equation (2) is modified by a 
diffusion term 

(Jf t ~D ^2 ) p ( w ' *) = K J du i P ( u i> *) / du ? p ( u 2, t) (18) 

x{S[v-u 1 + (l-p)g] -S(v-u 2 )} 

The temperature changes according to dT/dt + AT 3 / 2 = 2D, so the steady 
state temperature is T x = (2Z?/A) 2 / 3 . The relaxation toward the steady state 
is exponential, \T — X^l ~ exp(— const, x t). 

At the steady state, the Fourier transform F oc (k) = F(k, oo) satisfies 

{l + Dk 2 )F 00 (k)=F 00 (k-pk)F 00 (pk). (19) 

with D = -D/V^oo- Conservation of the total number of particles and the total 
momentum impose Foo(0) = 1 and F^O) = 0, respectively. The solution is 
found recursively to give the following infinite product [40] 

oo I _ l l \ 

^oo(fc) = n n [i + p 2m (i - v? (i - m) Dk 2 ] uj . (20) 

1=0 m=0 

Thus in one dimension the Fourier transform in determined analytically in the 
steady state. To extract the high-energy tail from (20) we note that F oa (k) 

has an infinite scries of poles located at ±i [p 2m (l - p) 2 ('- m ) D] 1 . The 
simple poles at ii/v^D closest to the origin imply an exponential decay of 
the velocity distribution [48,52], 

P.M-^e*', with v* = VD, (21) 

when \v\ — > oo. A re-summation similar to that used in [40] can be performed 
to evaluate the residue to this pole, and in turn, the prefactor 

2 When the mean-free path is comparable with the system size, the boundary 
effectively plays the role of a thermal heat bath. 
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This prefactor rapidly diverges, A(p) oc ^Jp exp [7T 2 /(12p)], as p — > 0. We 
have obtained only the leading asymptotic behavior. It will be interesting to 
determine for what range of velocities this behavior actually sets in. 

The leading high-energy behavior can also be derived by using a useful 
heuristic argument [26,28,48]. For sufficiently large velocities, the gain term 
in the collision integral in Eq. (18) is negligible. The resulting equation for 
the steady state distribution 

D^Poo(v) = -Poo(v) (23) 

yields the exponential high-energy tail (21). This argument applies to ar- 
bitrary collision rates. For example, if K oc v s , the right-hand side in (23) 
becomes —v s P oc implying that Poo{v) oc cxp(— |t>| 7 ) with 7=1 + 5/2. For 
hard spheres (5=1) one finds 7 = 3/2 [28], and curiously, the Gaussian tail 
arises only for the so-called very hard spheres (5 = 2) [35]. 

Finally, we notice that steady state properties in the heated case are 
intimately related with the relaxation properties in the cooling case. This 
can be seen via the cumulant expansion 



Foo(k) = exp 



(24) 



Replacing the term (1 + Dk 2 ) with exp[— n _1 (— Dfc 2 )"], and substi- 
tuting the cumulant expansion into Eq. (19) yields ip n = [na2n(p)}~ 1 ■ The 
cumulants defined via lnFoo(fc) = J2 n =i K n (ik) n /n\. Therefore, the 

steady-state cumulants arc directly related to the relaxation coefficients (12), 
= ^ 2 "~ 1 ' |! ^- (the odd cumulants vanish). 



3 Uniform Gases: Arbitrary Dimension 
3.1 The freely cooling case 

In general dimension, the colliding particles exchange momentum only along 
the impact direction. The post-collision velocities vi^ are given by a linear 
combination of the pre-collision velocities U1.2, 

vi,2 = ui,2 T (1 -p) (g • n)n. (25) 

Here g = Ui — U2 is the relative velocity and n the unit vector connecting the 
particles' centers. The normal component of the relative velocity is reduced 
by the restitution coefficient r = 1 — 2p, and the energy dissipation is given 
by AE = -p(l -p)(g • n) 2 . 

In random collision processes, both the impact direction and the iden- 
tity of the colliding particles arc chosen randomly. In such a process, the 
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Boltzmann equation 

= K J dn J duiPiui't) Jdu 2 P(u 2 ,t) (26) 
x{6 [v - ui + (1 - p)(g • n)n] - <5(v - m)} 

exactly describes the evolution of the velocity distribution function P(v,t). 
The overall collision rate is chosen to represent the typical relative velocity, 
K = \/T, with the granular temperature now being the average velocity 
fluctuation per degree of freedom, T = ^ J <ivu 2 P(v, t) with v = |v|. The 
evolution equation involves integration over all impact directions, and this 
angular integration should be normalized, J dn = 1. We tacitly ignored the 
restriction g • n > on the angular integration range in Eq. (26) since the 
integrand obeys the reflection symmetry n — > n. 

Several temporal characteristics such as the temperature and the autocor- 
relation behave as in the one-dimensional case. For example, the temperature 
satisfies Eq. (3) with prefactor A = 2p(l —p)/d reduced 3 by a factor d. The 
temperature therefore decays according to Haff's law (4), with the time scale 
to = dj [p(l — p) VTo] set by the initial temperature. Similarly, the decay rate 
of the autocorrelation function is merely reduced by the spatial dimension, 
-j^A(t' , t) = — -j 2 A(t',t). Consequently, the nonuniversal decay (16) and 
the logarithmic spread (17) hold in general. 

The Fourier transform, F(k, t) = J dv e lk v P(v, t), satisfies 

|- F(k, r) + F(k, r) = JdnF(k-p, r) F(p, r) (27) 

with p = (1 — p)(k • n) n reflecting the momentum transfer occurring during 
collisions. This equation was obtained by multiplying Eq. (26) by e tk v and 
integrating over the velocities. The power of the Fourier transform is even 
more remarkable in higher dimensions 4 as it reduces the (3d— 1)— fold integral 
in Eq. (26) to the (d - l)-fold integral in Eq. (27). 

Hereinafter, we consider only isotropic velocity distributions. The Fourier 
transform depends only on k = |k|, so we write F(k, r) = F(y, r) with y = k 2 . 
To perform the angular integration, we employ spherical coordinates with 
the polar axis parallel to k, so that k • n = cosd. The ^-dependent factor of 
the angular integration measure dm is proportional to (sin0) d ~ 2 d8. Denoting 

3 The reduction in rate by a factor d is intuitive because of the d independent 
directions only the impact direction is relevant in collisions. Mathematically, the 
prefactor A = 2p(l— p) J dnn\ is computed by using the identity n\ + . . - +n 2 d = 1 
that yields J dn n\ — 1 /d. 

4 For elastic Maxwell molecules, the Fourier transform was first used in unpublished 
theses by Krupp [34] , and then rediscovered and successfully utilized by Bobylev 
(see [35,36] for a review). 
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angular integration with brackets, (/) = J dnf, and using /j, = cos 2 9 gives 

(/)= I' d ^~lf~l\7" f(»), (28) 
Jo 

where B(a, b) is the beta function. This integration is properly normalized, 
(1) = 1. The governing equation (27) for the Fourier transform can now be 
rewritten in the compact from 

F(y, t) + F(y, r) = ( F(£y, r) F(rjy, r) ) , (29) 

with the shorthand notations £ = 1 — (1 — p 2 )[i and 77 = (1 — p) 2 ^. Unlike the 
one-dimensional case, explicit solutions of this nonlinear and nonlocal rate 
equation are cumbersome and practically useless. Nevertheless, most of the 
physically relevant features of the velocity distributions including the large 
velocity statistics and the time dependent behavior of the moments can be 
obtained analytically. 

We seek a scaling solution: P(v,t) = T~ d l 2 V(w) with w = vT" 1 / 2 , or 
equivalently F(y, r) = <P(x) with x = yT. The scaling function <P(x) satisfies 

-\x$'(x) + <P(x) = ${nx)) (30) 

and the boundary condition <P(x) = 1 — | x + ■ ■ ■ as x — > 0. In the elastic 
case, the velocity distribution is purely Maxwellian, <P(x) = er x l 2 . Indeed, 
A = and £ + r/ = 1 in this case. A stochastic process of elastic collisions 
effectively randomizes the velocities and leads to a thermal distribution [31]. 
In practice, this collision algorithm is used in Molecular dynamics simulations 
to thermalize velocities. 

From the one-dimensional case, we anticipate that the velocity distribu- 
tion has an algebraic large velocity tail. Generally, the large velocity behavior 
can be determined from the small wave number behavior of the Fourier trans- 
form. For example, the small- x expansion of the one-dimensional solution (10) 
contains both regular and singular terms: <£(x) = 1 — | x + ^ a; 3 / 2 + ■ ■ ■, and 
the dominant singular x 3 / 2 term reflects the w~ 4 tail of V(w). In general, if 
V(w) has an algebraic tail, 

V(w) ~ w~ a as w — > 00, (31) 

then <P(x) oc / °° dww d ^ 1 V(w) e lw ^ contains, apart from regular terms, the 
following dominant singular term: <£ s ing(a;) ~ a;(°' _d )/ 2 as x — > 0. The ex- 
ponent a can be now obtained by inserting <P(x) = ^ rcg (a;) + ^ s ing(a;) into 
Eq. (30) and balancing the dominant singular terms. We find that a is a root 
of the integral equation 5 

l-\1—^ = (^- d )/ 2 +^- d )/ 2 ). (32) 
5 This result was derived independently in [43,46]. 
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Fig. 1. The exponent a versus the collision parameter p for d — 2, and 3. For 
comparison with the leading large-dimension behavior (33) , a is rescaled by d. 



This equation can be recast as the eigenvalue problem A„ = v\\ using 
A„ = (1 — — rf) and v = (a — d)/2. It can also be expressed using the 
hypcrgeornetric function 2-Fi(d, b; c; z) [60] and the Euler's gamma function: 



l-p(l-p) 



a — d 
d 



2 -Pi 



d — a 1 d 
2 '2' 2 



(i-pr 



-d+l^ 



r 



Clearly, the exponent a = <r(d 7 p) depends in a nontrivial fashion on the 
spatial dimension d as well as the collision parameter p. 

There are two limiting cases where the velocity distribution approaches a 
Maxwellian, thereby implying a divergence of the exponent a. In the quasi- 
elastic limit (p — ► 0) we have a — > d/p. Even a minute degree of dissipation 
strongly changes the character of the system, and hence, energy dissipation 
is a singular perturbation [24,29]. As the dimension increases, the relative 
weight of the impact direction diminishes, and so does the role of dissipa- 
tion. For large dimensions, the integral {rf" - ^^) in Eq. (32) vanishes expo- 
nentially and asymptotic analysis of the remaining integral shows that the 
exponent grows linearly with the dimension: 



rf 



1 + \v - P 3 



,1/2 



.1/2 



p{l-p 2 ) 



(33) 



when d — > oo [43]. Equation (33) provides a decent approximation even at 
moderate dimensions (Fig. 1). Overall, the exponent <r(d,p) increases mono- 
tonically with increasing rf, and additionally, it increases monotonically with 



12 Ben-Nairn and Krapivsky 



decreasing p. Both features are intuitive as they mirror the monotonic de- 
pendence of the energy dissipation rate A = 2p{l — p)/d on d and p. Re- 
markably the exponent is very large. The completely inelastic case provides 
a lower bound for the exponent, a(d,p) > <r(d, 1/2) with a(d, 1/2) = 6.28753, 
8.32937, for d — 2, 3, respectively. For typical granular particles, a(d = 3,p = 
0.1) = 30, and such algebraic decays are impossible to measure in practice. 

The algebraic tail of the velocity distribution implies that moments of the 
scaling function <P(x) with sufficiently large indices diverge. In the scaling 
regime, moments of the velocity distribution can be calculated by expanding 
the Fourier transform in powers of x, 

M 

^SE^M" ( 34 ) 

71 = 

The order of terms in this expansion must be smaller than the order of the 
singular term, x v . The coefficients <p n , needed for calculating transport coef- 
ficients [53], yield the leading asymptotic behavior of the velocity moments, 
M k (t) = J dvv k P(v,t), via the relation (2n)\T n <f> n ~ {fi n )M 2n . Inserting 
the moment expansion into (30) yields a closed hierarchy of equations 

ri-1 

(A„ - n\l)<j) n = ^ ^m,n-m<Pm<l>n- m , (35) 
m—1 

with A„ = (1 — £™ — 77") and X m j = (£ m r/). Calculation of these coefficients 
requires the following integrals 

,„nv = r d) r ( n +h) = 1_L_ 2(n-l) + l 
V T(i)r(n+f) dd + 2 2(n-l) + d" 

Of course, (fio = 1 and <j>i = 1/2; further coefficients can be determined 
recursively from (35), e.g., 

1 l-3±r 



^ini' (36) 

^ ° d+2 

This coefficient is finite only when A2 > 2Ai or d < d^ = 1 + 3p 2 . Generally 
4> n is finite only when the left hand side of Eq. (35) is positive, A„ > nAi, or 
cquivalcntly, when the dimension is sufficiently small, d < d n . The crossover 
dimensions d n are determined from A„ = n\\. 

For d > d n , moments with index smaller than 2n are characterized by the 
temperature, while higher moments exhibit multiscaling asymptotic behavior. 
The time evolution of the moments can be studied using the expansion 

00 

F(y,r) = J2fn(r)(-y) n . (37) 

n=0 
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The actual moments are related to the /„ coefficients via (2n)!/„ = (fj, n }M 2n . 
Substituting the expansion (37) into (29) yields the evolution equations 

d ™ _1 

~^fn "f" A n / n ^ ^ ^m,n—mfmfn—m- (^8) 
m— 1 

We have /i oc e~ Air , this is just the HafT's law. From df 2 /dr + X 2 f 2 = M,ift 
we see that f 2 is a linear combination of two exponentials, e~ A2T and e~ 2AlT . 
The two decay coefficients are equal A2 = 2Ai at the crossover dimension 
d = d 2 . As expected, when d < d 2 the fourth moment is dictated by the 
second moment, f 2 oc f 2 . Otherwise, f 2 oc cxp(— A 2 r). In general, 

M 2n oc ( eXp ;~f j<j" ; (39) 
[ exp(-A„r) d> d n . 

Fixing the dimension and the collision parameter, moments of sufficiently 
high order exhibit multiscaling asymptotic behavior. In practice, the expo- 
nent a is typically large, and multiscaling occurs only for very high order 
moments. 



3.2 The forced case 



We consider white noise forcing as in the one-dimensional case. The steady 
state Fourier transform, F^lk) = F oc (k 2 ), satisfies 



(l + Dk 2 )F 0O (k 2 ) = (F 0O (^k 2 )F 0O (r 1 k 2 )). 



(40) 



This equation is solved recursively by employing the cumulant expansion 



(24). Writing 1 + Dk 2 = cxp ^2 n>1 {-Dk 2 ) n /n , we recast Eq. (40) into 



1 = 



exp 



n=l 



(41) 



with the auxiliary variables tp n = i/j n (l — £ n — f] n ). The coefficients ip n are 
obtained by evaluating recursively the angular integrals of the auxiliary vari- 
ables, (tpn), an d then using the identities tp n = (tp n )/X n . The few first coeffi- 
cients can be determined explicitly, e.g., V>i = 1/A and 



^2 



3d 2 



4(rf + 2)(l-p 2 )-12(l-p) 2 (l+p 2 )- 



(42) 



The nonvanishing second cumulant shows the steady state distribution is not 
Maxwellian. Moreover, the poles of the Fourier transform located at k = 
±i/^/D indicate that the large velocity tail of the distribution is exponential 
as in the one-dimensional case (21). Exponential decay is also suggested by 
the heuristic argument detailed above. 




Fig. 2. The correlation measure U for spatial dimension d — 3. Shown are the freely 
cooling case and the forced case. 



3.3 Velocity correlations 

Maxwellian velocity distributions were originally obtained for random elastic 
collision processes (see Ref. [32], p. 36). Maxwell's seminal derivation involves 
two basic assumptions: (1) The velocity distribution is isotropic, and (2) Cor- 
relations between the velocity components are absent. In inelastic gases, the 
velocity distributions arc non-Maxwcllian therefore, there must be correla- 
tions between the velocity components 6 . The quantity 

(vlvl)-{vl){v 2 v ) , 

provides a natural correlation measure. A non-vanishing U indicates that 
velocity correlations do exist, and the larger U the larger the correlation. 

9 
)k 

In the freely cooling case, U = 802 — 1; in the forced case, 



To compute U, we apply the identities: (v x ) = ^p- F 

»*2 »*S * k=0 

U = ip2/i>i- Substituting the corresponding coefficients yields 



and {v 2 x vl) = 

k=0 " 



6p 2 

d - (1 + 3p 2 ) ' 



C/cooiing = - — n : ^ 2N , (44a) 



f/forced " (d + 2)(1 + p) - 3(1- p)(l + p2) ' ( 44b ) 

6 Inelastic collisions discriminate the impact direction, thereby generating correla- 
tions among the velocity components. 
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Generally, correlations increase monotonically with p. Correlations are much 
weaker in the presence of an energy source because the nature of the driving 
is random (Fig. 2). The perfectly inelastic case provides an upper bound. For 
example, when d = 3 and p = 1/2 we have U coo u ng = 6/5 and [/forced = 2/15. 
Correlations decay as U oc p 2 and U k in the respective limiting cases 
of p — ► and d — > oo, where Maxwellian distributions are recovered. 

4 Impurities 

The dynamics of an impurity immersed in a uniform granular fluid often defies 
intuition. The impurity does not generally behave as a tracer particle, and 
instead it may move either faster or slower than the fluid. Despite extensive 
studies, many theoretical and experimental questions regarding the dynamics 
of impurities remain open [61-67]. 

Impurities represent the simplest form of polydispersity, an important 
characteristic of granular media [68,69]. Theoretically the impurity problem 
is a natural first step in the study of mixtures as it involves fewer collision 
parameters. The fluid background is not affected by the presence of the im- 
purity, and the impurity can be seen as "enslaved" to the fluid background. 

We study dynamics of a single impurity particle in a uniform background 
of identical inelastic particles. We set the fluid particle mass to unity and 
the impurity mass to m. Collisions between two fluid particles arc character- 
ized by the collision parameter p as in Eq. (25), while collisions between the 
impurity and any fluid particle are characterized by the collision parameter 
q. When an impurity particle of velocity Ui collides with a fluid particle of 
velocity u 2 its post-collision velocity vi is given by 



Two restitution coefficients, r p = r = 1 — 2p and r q = m — (m + l)q, char- 
acterize fluid-fluid and impurity-fluid collisions, respectively. The restitution 
coefficients obey < r < 1, and the collision parameters accordingly satisfy 
< p < 1 /2 and < q < -^-^ . The energy dissipated in an impurity- fluid 
collision is AE = -m(l -q)[2-(m+ 1)(1 - q)](g ■ n) 2 . 

Let Q(v, t) be the normalized velocity distribution of the impurity. In a 
random collision process, the impurity velocity distribution evolves according 
to the linear Lorentz-Boltzmann equation 



while the fluid distribution obeys (26). Two rates, K v = K and K q , char- 
acterize fluid-fluid and fluid-impurity collisions, respectively. These rates can 



Vi = Ui - (1 - q) (g • n) n. 



(45) 



dQ(y,t) 
dt 
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be eliminated from the respective equations (26) and (46) by introducing 
the collision counters, t p = t = J Q dt' K p (t'), and r q = J Q dt' K q (t'). The 
Lorentz-Boltzmann equation is further simplified by using the Fourier trans- 
form, G(k,t) = Jdve lkv Q(v,t) 

A G(k) + G(k) = J dnG(k - q) F(q), (47) 

with q = (1 — q) (k- n) n. This equation supplements the fluid equation (27). 

We consider two versions of the Maxwell model: An idealized case with 
equal collision rates, K p = K q (model A); and a more physical case with 
collision rates proportional to appropriate average relative velocities (model 
B). In the forced case, one can obtain the impurity velocity distribution using 
the cumulant expansion, and generally, velocity statistics of the impurity 
are similar to the background. Below, we discuss the more interesting freely 
cooling case. 



4.1 Model A 

When the two collision rates are equal, all pairs of particles are equally likely 
to collide with each other. For simplicity we set the overall rate to unity, 
K p = K q = 1; then both collision counters equal time, t p = r q = t. 

In the freely cooling case, the fluid temperature decays exponentially, 
T(t) = To exp [— 2p(l — p)t/d]. The impurity temperature, 0(t), defined by 
0(t) = ^ J dvv 2 Q(v, t), is coupled to the fluid temperature via a linear rate 
equation 

^e = - l —^e + ^-^T. (48) 

at d d 

The solution to this equation is a linear combination of two exponentials 

9(t) = (0 O - cT ) e - (1 - q ^^ d + cT e -Mi-p)t/d m ( 49) 

The constant c = (1 — q) 2 /[l — q 2 — 2p(l — p)] is simply the ratio between the 
impurity temperature and the fluid temperature in the long time limit. This 
generalizes the elastic fluid result c=(l — q)/(l + q) [70,71]. Generally, the 
impurity and the fluid have different energies, and this lack of equipartition 
is typical to granular particles [6,72] 

There are two different regimes of behavior. When 1 — q 2 > 2p(l —p), the 
impurity temperature is proportional to the fluid temperature asymptotically, 
— > c as t — > oo. In the complementary region 2p(l — p) > 1 — q 2 the 
ratio of the fluid temperature to the impurity temperature vanishes. Since the 
system is governed by three parameters (m, r p , r q ), it is convenient to consider 
the restitution coefficients as fixed and to vary the impurity mass. From the 
definition of the restitution coefficients the borderline case, q 2 = 1 — 2p(l— p), 
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defines a critical impurity mass 

r q + J(l + r 2 )/2 

m* = V , (50) 

1-^/(1 + ^/2 

This critical mass is always larger than unity. Asymptotically, the impurity 
to fluid temperature ratio exhibits two different behaviors 

gft) _^ f c m < m*, 

T(t) ~^ I oo m > m*. ^ ' 

We term these two regimes, the light impurity phase and the heavy impurity 
phase, respectively. In the light impurity phase, the initial impurity tem- 
perature becomes irrelevant asymptotically, and the impurity is governed by 
the fluid background. In the heavy impurity phase, the impurity is infinitely 
more energetic compared with the fluid and practically, it sees a static fluid. 
Interestingly, the dependence on the dimension is secondary as both the crit- 
ical mass, , and the temperature ratio c are independent of d. Below we 
study velocity statistics including the velocity distribution and its moments 
primarily in one dimension, where explicit solutions are possible. 



The light impurity phase In this phase, m < , the impurity is governed 
by the fluid background. We therefore seek scaling solutions of the same form 
as the fluid's: Q(v,t) = T~ l / 2 Q (uT -1 / 2 ) and G(k,t) = g(\k\T^ 2 ). From 
the Lorentz-Boltzmann equation (47), the latter scaling function satisfies the 
linear equation 

-p(l - p)zg'(z) + g(z) - g(qz) [1 + (1 - q)z] e~^> . (52) 

Since the fluid scaling function is a combination of z n e~ z with n — and 
n = 1, we try the series ansatz 

oo 

g(z) = Y,9nZ n e- z (53) 

71=0 

for the impurity. The first few coefficients, go = gi = 1 and g 2 = (1 — c)/2, 
follow from the small- z behavior g(z) = 1 — \cz 2 . The rest are obtained 
recursively 

q n - 1 (l-q)-p(l-p) 
9n = ~, T, v— ffn-1- 54 

1 — q n — np(l — p) 

The Fourier transform can be inverted to obtain the impurity velocity dis- 
tribution function explicitly. The inverse Fourier transform of e~ KZ is — ■> f ■> : 
the inverse transforms of z n e" z can be obtained using successive differen- 
tiation with respect to n. Therefore, the solution is a series of powers of 
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Lorcntzians 

2 °° 1 
2H%E^ (^) 

The coefficients Q n are linear combinations of the coefficients g^s with k < 
n + 1, e.g., Q 2 = 1 - 3.g2 + 3g 3 and Q3 = 4g 2 - 24g 3 + 6O54. There are 
special values of q for which the infinite sum terminates at a finite order. Of 
course, when q = p, the impurity is identical to the fluid and Q n = for all 
n > 2. When q 2 (l — q) = p(l — p), one has Q n = for all n > 3, so the 
velocity distribution includes only the first two terms. Regardless of q, the 
first squared Lorentzian term dominates the tail of the velocity distribution 
Q(w) ~ V(w) <~ w~ 4 , as w — > 00. One can show that this behavior extends 
to higher dimensions, Q(u>) ~ V(w) ~ w~ CT . Therefore, the impurity has the 
same algebraic extremal velocity statistics as the fluid. 

While the scaling functions underlying the impurity and the fluid are 
similar, more subtle features may differ. Moments of the impurity velocity 
distribution, L n (t) = j dvv n Q(v,t), obey the recursive equations 



, n-1 

-L n + b n L n =J2[ U )l m ( 1 - q) n - m L m M n - m , (56) 



dt 

m=2 



with b n (q) = 1 — q n . Asymptotically, the fluid moments decay exponentially 
according to M n {t) cx e ^ a "(p) t . Using this asymptotic behavior, we ana- 
lyze the (even) impurity moments. The second moment, i.e. the impurity 
temperature, was shown to behave similar to the fluid temperature when 
a 2{p) < b2(q). The fourth moment behaves similarly to the fourth moment 
of the fluid, L4 cx M4, when a 4 (p) < 64(9). However, in the complementary 
case, the fourth moment behaves differently, L n (t) cx e~ b4( ^*. In general, 
when a n (p) < b n (q), the nth impurity moment is proportional to the nth 
fluid moment, L n cx M n . Otherwise, when a n {p) > b n (q), the nth impurity 
moment is no longer governed by the fluid and 

L n (t) cx e^ 1 -^. (57) 

This series of transition affecting moments of decreasing order occurs at 
increasing impurity masses, 

mi > m 2 > ■ • ■ > "ioo- (58) 

When m > m n , the ratio M 2 k/L 2 k diverges asymptotically for all k > n. The 
transition masses 

m n = ^ + H(l-^) 2 " + (l + ^) 2 1- (59) 
1 - - r P ) 2n + (1 + r p ) 2n }^ 

are found from q 2n = p 2n + (1 — p) 2n and the definitions of the restitu- 
tion coefficients. All of the transition masses are larger than unity, so the 
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impurity must be heavier than the fluid for any transition to occur. The 
largest transition mass is mi = m*, and the smallest transition mass is 
moo = lim„^oo m n = (1 + r. p + 2r q )/(l — r p ). Impurities lighter than the 
latter mass, m < moo, mimic the fluid completely 

The heavy impurity phase In this phase, m > m*, the velocities of the 
fluid particles are asymptotically negligible compared with the velocity of the 
impurity. Fluid particles become stationary as viewed by the impurity, and 
effectively, u 2 = in the collision rule (45): 



Mathematically this process is reminiscent of a Lorentz gas [73]. Physically, 
the two processes are different. In the granular impurity system, a heavy par- 
ticle scatters of a static background of lighter particles, while in the Lorentz 
gas the scatterers are infinitely massive. 

We first consider the one-dimensional case. Setting ui = in the Lorcntz- 
Boltzmann equation (46), integration over the fluid velocity 112 is trivial, 
/ du2 P{u2, t) = 1, and integration over the impurity velocity u\ gives 



This equation can be solved directly by considering the stochastic process the 
impurity particle experiences. In a sequence of collisions, the impurity velocity 
changes according to v — > qv — ► q 2 v — > ■ ■ ■ with w the initial velocity. 
After n collisions the impurity velocity decreases exponentially, v„ — q n vo. 
Furthermore, the collision process is random, and therefore, the probability 
that the impurity undergoes exactly n collisions up to time t is Poissonian 
t n e~ l jn\. Thus, the velocity distribution function reads 



where Qo(v) is the initial velocity distribution of the impurity. The impurity 
velocity distribution function is a time-dependent combination of "replicas" 
of the initial velocity distribution. Since the arguments are stretched, compact 
velocity distributions display an infinite set of singularities, a generic feature 
of the Maxwell model. 

In contrast to the velocity distribution, the moments L n (t) exhibit a much 
simpler behavior. Indeed, from Eq. (61) one finds that every moment is cou- 
pled only to itself, ^L n = —(1 — q n )L n . Solving this equation we recover 
Eq. (57); in the heavy impurity phase, however, it holds for all n. There- 
fore the moments exhibit multiscaling asymptotic behavior. The decay coef- 
ficients, characterizing the n-th moment, depend on n in a nonlinear fashion. 



v = u — (1 — q) (u • n) n. 



(60) 




(61) 




(62) 
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4.2 Model B 



For hard sphere particles, the collision rates are proportional to the relative 
velocity. The overall collision rates in the Maxwell model represent the aver- 
age relative velocity, and a natural choice is ((vi — v 2 ) 2 ) oc {T\ + T 2 )/2. 
Therefore, the rates K p = \ff and K q = y/{T + &)/2, should be used in the 
Boltzmann equation (26) and the Lorentz-Boltzmann equation (46), respec- 
tively. This modification suppresses the heavy impurity phase, although the 
secondary transitions corresponding to higher order moments remain. 
The fluid and the impurity temperatures obey 



dt 



-T 



-VT 



2p(l-p) 



T 



±e = J T+0 

dt 



d 



9 + 



(1-9) 2 



Consequently, the temperature ratio, S = O/T, evolves according to 

-.2 



VT dt 



1 + S 



1 



d 



d 



+ Ml - P) S. (63) 
d 



The loss term, which grows as S 3 ^ 2 , eventually overtakes the gain term that 
grows only linearly with S. Therefore, the asymptotic temperature ratio re- 
mains finite, S — > c, where c is the root of the cubic equation 



1 + c 



1-g 
l + Q 



2p(l-p ) 
1-q 2 



■ C. 



(64) 



Consequently, there is only one phase, the light impurity phase. Intuitively, 
as the impurity becomes more energetic, it collides more often with fluid 
particles. This mechanism limits the growth rate of the temperature ratio. As 
in model A, energy equipartition does not generally occur, and the behavior 
is largely independent of the spatial dimension. 

Qualitatively, results obtained for the light impurity phase in model A ex- 
tend to model B. The impurity velocity distribution follows a scaling asymp- 
totic behavior. The only difference is that the collision terms are proportional 
to = + c)/2, and Eq. (52) generalizes as follows 



-(3p{\ - p)zg'(z) + g(z) = (1 + az) e~ az g{qz). 



(65) 



Seeking a series solution of the form (53) leads to the following recursion 
relations for the coefficients 



g n = 



( 1 - q)q n - 1 - pp(l - p) 
1 — q n — nf3p(l — p) 



(66) 



Again, the velocity distribution is a combination of powers of Lorentzians 
as in Eq. (55). Moreover, the coefficients Q n are linear combinations of the 
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coefficients g„'s as in model A. Most importantly, the large- velocity tail is 
generic Q(w) ~ w~ 4 . 

The fluid moments evolve according to (11) with r = t p , while the impu- 
rity moments evolve according to 



Note that asymptotically r q — > t p //3 7 so the fluid moments decay according 
to M n cx e ~' 3a «(p) T 'j. Hence, when (3a n (p) < b n (q), the impurity moments 
are enslaved to the fluid moments, i.e., L n cx M n asymptotically. Other- 
wise, sufficiently large impurity moments behave differently than the fluid 
moments, viz. M n cx e~ b " T «. Although the primary transition affecting the 
second moment docs not occur (mi = to* = oo), secondary transitions af- 
fecting larger moments do occur at a scries of masses, as in Eq. (58). The 
transition masses m„ are found by solving /3a n (p) = b n (q) simultaneously 
with Eq. (64). For example, for completely inelastic collisions (r p = r q = 0) 
one finds m 2 = 1.65. These transitions imply that some velocity statistics of 
the impurity, specifically large moments, are no longer governed by the fluid. 

4.3 Velocity Autocorrelations 

The impurity autocorrelation function satisfies dA/dr q = —(1 — q)A. There- 
fore, its decay is similar to (15) 



with q replacing p and r q replacing r = r p . For model A where the collision 
counters equal time, the decay remains exponential. However, the impurity 
autocorrelation decays with a different rate than the fluid. For model B, one 
can show that the algebraic decay (16) holds asymptotically with p replaced 
by q. We conclude that while one-point velocity statistics of the impurity are 
governed by the fluid, two-time statistics are different. 

5 Mixtures 

Granular media typically consists of mixtures of granular particles of several 
types. In contrast with the impurity problem, the different components of a 
mixtures are coupled to each other. Moreover, there are numerous param- 
eters, making mixtures much harder to treat analytically. Remarkably, the 
impurity solution can be generalized to arbitrary mixtures 7 . We detail the 
freely cooling case in one dimension. 

7 Mixtures were treated analytically in the elastic Maxwell model [74] 




(67) 



A(T q ,T q ) - A,CX P [-(1 - q){ Tq - Tg)] 



(68) 
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Consider a binary mixture where particles of type 1 have mass ra\ and 
concentration c\, and similarly for particles of type 2. We set C\ + c 2 = 1. 
Also, unit collision rates are considered for simplicity. Collisions between a 
particle of type i = 1, 2 and a particle of type j = 1,2 are characterized by 
the collision parameter p^ — (rrii — Tij ) / (m t + mj ) , with the restitution 
coefficient. We denote the normalized velocity distribution of component i by 
Pi(v,t) and its Fourier transform by Fi(k,t). The governing equations now 
couple the two distributions 

d 2 

— F i (k)+F i (k)=^2c j F j (k- Pij k) Fi { Pij k) . (69) 
i=i 

Let Ti = J dv v 2 Pi(v, t) be the temperature of the ith component. Writing 
Fi(k, t) = 1 — ^k 2 Ti, the temperatures evolve according 



±(Ti\ = _ f An A 12 \ /Tj 
dt\T 2 { A 2 i A 22 T 2 



(70) 



The diagonal matrix elements are Xu = 2ciPu(l — pu) + Cj(l — p? ) with 
j 7^ i, and the off-diagonal matrix elements are Ajj = — Cj (1 —p^) 2 . Therefore, 
the temperatures are sums of two exponential terms, exp(— A±t). The decay 
coefficients are the two (positive) eigenvalues 

1 



A± = 2 



An + A 22 ± V( A n ~ A 2 2 ) 2 + 4Ai 2 A 2 i J . (71) 
Asymptotically, the term with the smaller decay rate A = A_ dominates 

IK^Ce-* with Q= (A + -W)-^(0), (?2) 

A_|_ — A_ 

The ratio between the two temperatures approaches C 2 /Ci = (A — An)/Ai2, 
and the two components have different temperatures [50]. Furthermore, as 
long as the components are coupled, both temperatures are finite. For ex- 
ample, a vanishing C 2 implies that one of the concentrations vanishes as 
A = A M . 

As in the cases of homogeneous gases and impurities, we seek a scaling 
solution of the form Pi(v,t) = e xt l 2 Vi{w) with w — ve xt / 2 , or equivalently, 
Fi(k,t) = fi(z) with z = \k\ e~ xt l 2 . The scaling functions fi(z) are coupled 
via the non-local differential equations 

j 2 

- g Xz fi( Z ) + fi ( Z ) = XI C 3 fj( Z ~ Pv Z ) fi (Pij Z )' ( 73 ) 

Substituting the series solution (53) , fi(z) — ^2^ =0 Ai^ n z n e~ z , yields recur- 
sion relations for the coefficients Ai^ n 

(\\ \ 2 n 

1 - — j Ai <n + -^ >n _i =Y,T, CjAj,mAi,n-mO--Pij) m Pv~ m - ( 74 ) 
' j=l m=0 
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The small z behavior fi(z) = 1 — \CiZ 2 implies the first three coefficients 
Ai t o = Ai : i = 1, and Ai^ = (1 — Cj)/2. For n > 3, the coefficients (Ai, n , A2, n ) 
are solved recursively in pairs. Each such pair satisfies two inhomogeneous 
linear equations. Therefore, as in the impurity case, the velocity distribu- 
tion is an infinite series of powers of Lorentzians (55). Although the two 
velocity distributions are different, they have the same extremal behavior, 
Vi(w) ~ u>~ 4 . It is straightforward to generalize the above to mixtures with 
arbitrary number of components, and to incorporate different collision rates, 
in particular, = y/(Ti + Tj)/2. 

6 Lattice Gases 

Inelastic collisions generate spatial correlations and consequently, inelastic 
gases exhibit spatial structures such as shocks. Thus far, we considered only 
mean-field collision processes where there is no underlying spatial structure. 
Random collision processes can be naturally generalized by placing particles 
on lattice sites and allowing only nearest neighbors to collide. 

Consider a one-dimensional lattice where each site is occupied by a single 
particle. Let Vj be the velocity of the particle at site j. The velocity of such a 
particle changes according to Eq. (1) due to interactions with either of its two 
neighbors. Time is conveniently characterized by the collision counter r. In 
an infinitesimal time interval dr, the velocity of a particle changes as follows 



This process is stochastic and we are interested in averages over all possible 
realizations of the process, denoted by an overline. We consider random ini- 
tial conditions where the average velocity vanishes and no correlations are 
present: Vj(0) = and Vi(0)vj(0) = T n 5ij. 

Spatial velocity correlations satisfy closed equations as in the Ising-Glauber 
spin model [75]. For example, from the dynamical rules (75), the average ve- 
locity Vj (r) = Vj (t) , obeys a discrete diffusion equation 



From the initial conditions, Vj(0) = 0, we obtain Vj(t) = 0. Consider the 
spatial correlation function vTwT. The initial state is translationally invariant, 
so this property persists. The correlation functions R n = VjVj+ n satisfy 




prob. I — 2dr: 
Vj-i(r)] prob. dr; 
v j+i( T )] prob. dr. 



(75) 



-± = (i-p)(y j - 1 -2V j + V j+1 ). 



(76) 



dR, 



2(1 -p)(R, 



2R n + Rn+l) 



n>2; 



dr 
dRi 

~dV 

dR 



'71—1 



2{l-p)\pR -(l+p)R 1 +R 2 }; 



(77) 
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The initial conditions are i?„(0) = T <5 nj o- Since we are interested in the 
asymptotic behavior, we employ the continuum approximation. The correla- 
tion function satisfies the diffusion equation dR/dr = 2(1 — p)d 2 R/d 2 n, so 
the solution is the Gaussian 



Rn[T) " VSd-pfrT 



n 2 



8(1 -p) T 



(78) 



The temperature, T = R decays as T(r) ~ T [8(1 — p)7tt] -1 / 2 in the long 
time limit 8 . 

Although no correlations were present in the initial conditions, spatial 
correlations develop at later times. The corresponding correlation length £ 
grows diffusively with the collision counter, £ <~ t 1 / 2 . The system consists of 
a network of domains of typical size £. Inside a domain, velocities are strongly 
correlated, and momentum conservation yields a relation between the average 
velocity and the domain size, ~ £ -1 / 2 ^ r -1 / 4 . This scale is consistent 
with the temperature behavior above, T ~ v 2 . In arbitrary dimension, this 
scaling argument yields T <~ T~ d l 2 [22,42]. At least for scalar velocities, the 
correlation functions obey closed equations in arbitrary dimension [76], and 
this behavior can be also obtained analytically. 

The actual time dependence is determined from the collision rate. We con- 
sider two choices. In model A, the rate is proportional to the typical velocity 
K 2 = T. In the physically more realistic model B, the rate is proportional 
to the average relative velocity K 2 = (vj+i — Vj) 2 = 2(Rq — R\) oc —dT/dr. 
Hence K - t~ q with a = d/4 (model A) and a = (d + 2)/4 (model B). The 
time t — J dr' / K{t') grows as t ~ r 1+Q , and therefore, 



-^-r model A; 

2d 
d+6 



*" 7 With 1=<W model B : ^ 



In either case, Haff 's cooling law is recovered only in the infinite dimension 
limit, d — > oo. Otherwise, the appearance of spatial correlations slows down 
the temperature decay. As neighboring particles become correlated their rel- 
ative velocity and consequently collision rate is reduced. Qualitatively similar 
behavior was shown for freely cooling inelastic gases, namely, breakdown of 
the mean-field cooling law due to the formation of strong spatial correlations 
among particles velocities [24,25]. 



7 Conclusions 

We presented analytic results for random inelastic collision processes. In gen- 
eral, when the collision rate is uniform, the convolution structure of the col- 
lision integrals translates to products in Fourier space. While the governing 

8 An exact solution of the discrete equations (77) is possible. It yields an identical 
asymptotic expression for the temperature. 
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equations are both nonlinear and nonlocal, they are closed and amenable to 
analytical treatment. In the freely cooling case, a small wave number analysis 
of the Fourier transform displays both regular and singular terms. The regu- 
lar terms yield the low order moments, while the leading singular term gives 
the high-energy tail. In the forced case, the Fourier transform is an analytic 
function and complex residue analysis yields the high-energy tail. 

In the freely cooling case, the velocity distribution approaches a scaling 
form and displays an algebraic large velocity tail. In one dimension, the scal- 
ing function is a universal; otherwise, it depends on the degree of inelasticity. 
The exponent governing the high-energy tail is typically very large, and may 
be difficult to measure in practice. We have also shown that sufficiently large 
moments exhibit multiscaling behavior, and hence, are not characterized by 
the temperature. The autocorrelation function decays algebraically, with an 
exponent that depends on the collision parameters. The spread of a tagged 
particle exhibits a universal Vint growth. 

We demonstrated that an impurity immersed in a uniform fluid may or 
may not mimic the background. Fixing the collision parameters, for suffi- 
ciently low impurity mass, the impurity's velocity distribution, velocity mo- 
ments and extremal velocity statistics are all governed by the fluid. However, 
there is a series of phase transitions occurring at a series of increasing masses, 
where impurity moments of decreasing orders decouple from the fluid. These 
transitions indicate that sufficiently heavy impurities are very energetic and 
effectively, they experience a static fluid background. 

For binary mixtures, we examined only primary velocity statistics and 
showed that qualitatively, all components have the same temperature decay 
and extremal velocity statistics. It remains to be seen whether the different 
components may exhibit different asymptotic behaviors of more subtle veloc- 
ity statistics such as the high order moments. The impurity case demonstrates 
how in certain limiting cases, mixtures may display anomalous behavior. 

In the forced case, injection of energy counters the energy dissipation and 
the system relaxes toward a steady state. We considered white noise forcing, 
and showed that the steady state distribution has an exponential high-energy 
tail. Steady state characteristics in the forced case arc directly related to time 
dependent relaxation characteristics in the freely cooling case. 

Results obtained in the framework of the inelastic Maxwell model are 
exact for random collision processes where spatial structure is absent, and 
particles collide irrespective of their relative velocities. Such results are an un- 
controlled approximation of the inelastic hard sphere problem. Nevertheless, 
as a conceptual tool, the inelastic Maxwell model is powerful. For example, 
it demonstrates the development of correlations among the velocity compo- 
nents, as well as spatial correlations. It also raises doubts concerning the 
suitability of several widely used techniques such as perturbation expansions 
in the quasi-elastic limit, and expansions in terms of Sonine polynomials. 
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We demonstrated that a lattice gas generalization of the Maxwell model 
remains tractable when the kinematic constraint (particles should collide only 
if they are moving toward each other) is ignored. Taking this constraint into 
account leads to shock-like structures in one dimension and to vortices in two 
dimensions [42]. Lattice models can therefore be used to study development 
of spatial correlations and spatial structures in inelastic gases. 

We presented in detail the most basic realization of the Maxwell model. 
Several other generalizations are feasible. Experiments can now measure the 
distributions of impact angles and of the effective restitution coefficients. Such 
phenomenological information can be incorporated into the rate equations. 
The integration measure can be redefined to give different weights to different 
angles. Moreover, a distribution of restitution coefficients can be introduced 
by integrating the collision integrals with respect to the collision parameters. 
For example, random collision processes were successfully used to model the 
role of the boundary in driven gases [77] . 

We restricted our attention to the case where all moments of the initial ve- 
locity distribution are finite. However, there are other infinite energy solutions 
of the Boltzmann equation. For example, any Lorentzian, F(z) = exp(— Cz), 
is a steady state solution of Eq. (6). A remaining challenge is to classify the 
evolution of an arbitrary velocity distribution. Also, it will be interesting to 
characterize the full spectrum of extremal behaviors, characterizing velocities 
far larger than the typical velocity. 

Finally, to model transport and other hydrodynamics problems one has 
to incorporate the spatial dependence explicitly. In the framework of the 
Maxwell model the simplest such problem - the shear flow - has been re- 
cently investigated by Cercignani [78]. For inelastic hard spheres, a number 
of unidirectional hydrodynamic flows were analyzed by Goldshtein and co- 
workers [79-81]. It would be interesting to investigate these problems in the 
framework of the Maxwell model. 
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